Session lead: Katie Smith (Centre for Ecology & Hydrology) k.a.smith@ceh.ac.uk
This is an demonstration of how to use the airGR package of hydrological models in R, as well as how to plot interactive timeseries graphs with the dygraphs package.
library(xts)
package <U+393C><U+3E31>xts<U+393C><U+3E32> was built under R version 3.3.3Loading required package: zoo
Attaching package: <U+393C><U+3E31>zoo<U+393C><U+3E32>
The following objects are masked from <U+393C><U+3E31>package:base<U+393C><U+3E32>:
as.Date, as.Date.numeric
library(dygraphs)
package <U+393C><U+3E31>dygraphs<U+393C><U+3E32> was built under R version 3.3.3
observed_data <- read.csv("~/Easter_Work/GR_teaching/Qobs_390010.csv")
observed_data
observed_data$DATE <- strptime(observed_data$DATE, format = "%d/%m/%Y")
observed_data_xts <- as.xts(observed_data$Qobs, order.by = observed_data$DATE)
dygraph(observed_data_xts)%>%
dyOptions()%>%
dyRangeSelector()
precip_data <- read.csv("~/Easter_Work/GR_Teaching/rain_1961_2014_390010.csv")
precip_data
PET_data <- read.csv("~/Easter_Work/GR_Teaching/CHESS_PET_1961_2015_390010.csv")
PET_data
precip_data$DATE <- strptime(precip_data$DATE, "%Y-%m-%d")
PET_data$DATE <- strptime(PET_data$DATE, "%Y-%m-%d")
first_date <- max(observed_data$DATE[1], precip_data$DATE[1], PET_data$DATE[1])
last_date <- min(observed_data$DATE[length(observed_data$DATE)], precip_data$DATE[length(precip_data$DATE)], PET_data$DATE[length(PET_data$DATE)])
# make an empty data frame
thames_data <- as.data.frame(matrix(NA,nrow=as.numeric((last_date-first_date)+1), ncol=4))
colnames(thames_data) <-c ("date","PET","precip","obs")
# make the date timeseries
thames_data$date <- seq(first_date, last_date, by="days")
# populate the data frame with the data
thames_data$obs <- observed_data$Qobs[which(observed_data$DATE==thames_data$date[1]):which(observed_data$DATE==thames_data$date[length(thames_data$date)])]
thames_data$precip <- precip_data$Mean_rainfall[which(precip_data$DATE==thames_data$date[1]):which(precip_data$DATE==thames_data$date[length(thames_data$date)])]
thames_data$PET <- PET_data$PET[which(PET_data$DATE==thames_data$date[1]):which(PET_data$DATE==thames_data$date[length(thames_data$date)])]
# convert the observed discharge to runoff (so its in the same units as the precip)
# divide by catchment area (m2) and mulitply by 86.4
thames_data$obs <- (thames_data$obs/9948.0)*86.4
thames_data_xts <- as.xts(thames_data[,3:4], order.by=thames_data$date)
# initiate the dygraph
dygraph(thames_data_xts, main = "Naturalised Runoff and Precipitation Observations for the Thames at Kingston")%>%
# define the first axis
dyAxis(dygraph = graphOut, name = "y", label = "runoff (mm/day)",
valueRange = range(thames_data_xts[, "obs"],
na.rm = TRUE)* c(0.01, 1.59))%>%
# define the second axis
dyAxis(dygraph = graphOut, name = "y2", label = "precip (mm/day)",
valueRange = rev(range(thames_data_xts[, "precip"],
na.rm = TRUE)* c(0.01, 2.99)))%>%
# plot the data
dySeries("obs",axis = 'y')%>%
dySeries("precip", axis = 'y2', stepPlot = TRUE,
fillGraph = TRUE)%>%
dyOptions(colors = RColorBrewer::brewer.pal(3,"Set1")[3:1]) %>%
dyRangeSelector()
require(airGR, quietly=TRUE)
package <U+393C><U+3E31>airGR<U+393C><U+3E32> was built under R version 3.3.3
BasinObs <- thames_data
colnames(BasinObs) <- c('DatesR','E','P', 'Qobs')
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J,DatesR = BasinObs$DatesR,
Precip = BasinObs$P,PotEvap = BasinObs$E)
str(InputsModel)
List of 3
$ DatesR : POSIXlt[1:19723], format: "1961-01-01 00:00:00" "1961-01-02 00:00:00" ...
$ Precip : num [1:19723] 8.6376 4.571 2.1015 0.0201 7.1859 ...
$ PotEvap: num [1:19723] 0.367 0.353 0.412 0.32 0.612 ...
- attr(*, "class")= chr [1:3] "InputsModel" "daily" "GR"
# note NA values of precip and PET are NOT ALLOWED
# first define the period to run the model over
Ind_Run <- seq(which(BasinObs$DatesR=="1981-01-01"),
which(BasinObs$DatesR=="2014-12-31"))
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
IndPeriod_WarmUp = NULL)
Model warm-up period not defined -> default configuration used
The year preceding the run period is used
Model states initialisation not defined -> default configuration used
str(RunOptions)
List of 6
$ IndPeriod_WarmUp: int [1:365] 6941 6942 6943 6944 6945 6946 6947 6948 6949 6950 ...
$ IndPeriod_Run : int [1:12418] 7306 7307 7308 7309 7310 7311 7312 7313 7314 7315 ...
$ IniStates : num [1:67] 0 0 0 0 0 0 0 0 0 0 ...
$ IniResLevels : num [1:2] 0.3 0.5
$ Outputs_Cal : chr "Qsim"
$ Outputs_Sim : chr [1:16] "DatesR" "PotEvap" "Precip" "Prod" ...
- attr(*, "class")= chr [1:3] "RunOptions" "GR" "daily"
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE,
InputsModel = InputsModel,
RunOptions = RunOptions,
Qobs = BasinObs$Qobs[Ind_Run])
str(InputsCrit)
List of 5
$ BoolCrit : logi [1:12418] TRUE TRUE TRUE TRUE TRUE TRUE ...
$ Qobs : num [1:12418] 0.635 0.623 0.587 0.571 0.566 ...
$ transfo : chr ""
$ Ind_zeroes: NULL
$ epsilon : NULL
- attr(*, "class")= chr "InputsCrit"
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel)
str(CalibOptions)
List of 3
$ FixedParam : logi [1:4] NA NA NA NA
$ SearchRanges : num [1:2, 1:4] 4.59e-05 2.18e+04 -1.09e+04 1.09e+04 4.59e-05 ...
$ StartParamDistrib: num [1:3, 1:4] 169.017 247.151 432.681 -2.376 -0.649 ...
- attr(*, "class")= chr [1:3] "CalibOptions" "GR4J" "HBAN"
OutputsCalib <- Calibration_Michel(InputsModel = InputsModel,
RunOptions = RunOptions,
InputsCrit = InputsCrit,
CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J,
FUN_CRIT = ErrorCrit_NSE)
Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
Screening completed (81 runs)
Param = 432.681 , -0.649 , 83.096 , 2.384
Crit NSE[Q] = 0.8923
Steepest-descent local search in progress
Calibration completed (19 iterations, 216 runs)
Param = 607.894 , -0.734 , 87.357 , 2.315
Crit NSE[Q] = 0.9246
Param <- OutputsCalib$ParamFinalR
Param
[1] 607.8936811 -0.7336304 87.3567230 2.3153153
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions,
Param= Param)
str(OutputsModel)
List of 16
$ DatesR : POSIXlt[1:12418], format: "1981-01-01 00:00:00" "1981-01-02 00:00:00" ...
$ PotEvap : num [1:12418] 1.347 0.38 0.795 0.57 0.454 ...
$ Precip : num [1:12418] 0.0972 0.1211 0.0523 0.1147 0.3248 ...
$ Prod : num [1:12418] 352 351 350 350 349 ...
$ AE : num [1:12418] 1.127 0.334 0.663 0.488 0.43 ...
$ Perc : num [1:12418] 0.387 0.383 0.378 0.374 0.372 ...
$ PR : num [1:12418] 0.387 0.383 0.378 0.374 0.372 ...
$ Q9 : num [1:12418] 0.355 0.35 0.345 0.341 0.338 ...
$ Q1 : num [1:12418] 0.0398 0.0393 0.0387 0.0382 0.0378 ...
$ Rout : num [1:12418] 42.2 42 41.7 41.4 41.2 ...
$ Exch : num [1:12418] -0.0592 -0.0577 -0.0564 -0.0551 -0.0539 ...
$ AExch : num [1:12418] -0.099 -0.097 -0.0951 -0.0933 -0.0917 ...
$ QR : num [1:12418] 0.599 0.578 0.559 0.542 0.526 ...
$ QD : num [1:12418] 0 0 0 0 0 ...
$ Qsim : num [1:12418] 0.599 0.578 0.559 0.542 0.526 ...
$ StateEnd: num [1:67] 374.4 45.5 NA NA NA ...
- attr(*, "class")= chr [1:3] "OutputsModel" "daily" "GR"
plot(OutputsModel, Qobs=BasinObs$Qobs[Ind_Run])
# make a few changes to the calibration criteria
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE2,
InputsModel = InputsModel,
RunOptions = RunOptions,
Qobs = BasinObs$Qobs[Ind_Run])
# rerun the calibration
OutputsCalib <- Calibration_Michel(InputsModel = InputsModel,
RunOptions = RunOptions,
InputsCrit = InputsCrit,
CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J,
FUN_CRIT = ErrorCrit_KGE2)
Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
Screening completed (81 runs)
Param = 432.681 , -0.649 , 83.096 , 2.384
Crit KGE'[Q] = 0.8589
Steepest-descent local search in progress
Calibration completed (51 iterations, 502 runs)
Param = 598.748 , -0.690 , 93.679 , 2.297
Crit KGE'[Q] = 0.9621
# redefine the parameters
Param <- OutputsCalib$ParamFinalR
# rerun the model
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions,
Param= Param)
# plot again
plot(OutputsModel, Qobs=BasinObs$Qobs[Ind_Run])
# go back to the beginning, redefine the period to run on (the period we haven't used for calibration, minus the first year needed for warm up)
Ind_Run <- seq(which(BasinObs$DatesR=="1962-01-01"),
which(BasinObs$DatesR=="1980-12-31"))
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
IndPeriod_WarmUp = NULL)
Model warm-up period not defined -> default configuration used
The year preceding the run period is used
Model states initialisation not defined -> default configuration used
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE2,
InputsModel = InputsModel,
RunOptions = RunOptions,
Qobs = BasinObs$Qobs[Ind_Run])
Param <- OutputsCalib$ParamFinalR
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions,
Param= Param)
OutputsCrit <- ErrorCrit_KGE2(InputsCrit = InputsCrit,
OutputsModel = OutputsModel)
Crit. KGE'[Q] = 0.9039
SubCrit. KGE'[Q] cor(sim, obs, "pearson") = 0.9466
SubCrit. KGE'[Q] sd(sim)/sd(obs) = 0.9253
SubCrit. KGE'[Q] mean(sim)/mean(obs) = 1.0282
Ind_Run <- seq(which(BasinObs$DatesR=="1962-01-01"),
which(BasinObs$DatesR=="2014-12-31"))
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
IndPeriod_WarmUp = NULL)
Model warm-up period not defined -> default configuration used
The year preceding the run period is used
Model states initialisation not defined -> default configuration used
Param <- OutputsCalib$ParamFinalR
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions,
Param= Param)
plot_output_data <- as.data.frame(matrix(NA, ncol = 3,
nrow = length(OutputsModel$DatesR)))
colnames(plot_output_data) <- c("Date", "Qsim", "Qobs")
plot_output_data$Date <- OutputsModel$DatesR
plot_output_data$Qsim <- OutputsModel$Qsim
plot_output_data$Qobs <- BasinObs$Qobs[Ind_Run]
plot_output_data_xts <- as.xts(plot_output_data, order.by = plot_output_data$Date)
library(RColorBrewer)
dygraph(plot_output_data_xts, main="Observed and Simulated Runoff for the Thames at Kingston (Naturalised)")%>%
dyOptions(colors = RColorBrewer::brewer.pal(3,"Set1")[3:1])%>%
dyAxis("y", label="Runoff (mm/day)")%>%
dyRangeSelector()